### Replicates Figure 1(A) for Study #2 in the paper ###

### Analysis ####
library(RItools)
library(xtable)
library(car)
setwd("~/Dropbox/Data/Paper-GDN/Analysis-Paper-01/Replic") #change to local folder
savedir <- getwd()
load("S02-dataBrazil.RData")

### Examine Randomization
table(d$TO,d$Tcond)

### Main effects of otherness
lm(OUTsupport~TO,data=d) 
lm(OUTsupport~TO,data=subset(d,attention==T)) 

### Main effect of conditionality
lm(OUTsupport~Tcond,data=d) #expected positive increase 
lm(OUTsupport~Tcond,data=subset(d,attention==T)) 

### Interaction effects
summary(lm(OUTsupport~Tcond*TO,data=d))
summary(lm(OUTsupport~Tcond*TO+attention,data=d))
summary(lm(OUTsupport~Tcond*TO,data=subset(d,attention==T)))


##### Figure 1(A) in the Paper ###############################

### A function for computing robsut SE ####
robse <- function(x,sig=0.95){
	require(sandwich)
	robSE <- sqrt(diag(vcovHC(x,type="HC2")))
	terms <- length(coef(x))
	robP <- rep(NA,terms)
	for(jj in 1:terms){
			robP[jj] <- linearHypothesis(x,diag(terms)[jj,],vcov=vcovHC(x,type="HC2"))$Pr[2]}
	upr <- coef(x) + qnorm(sig) * robSE
	lwr <- coef(x) - qnorm(sig) * robSE
	out <- cbind(robse=robSE,upr=upr,lwr=lwr,robpval=robP)
	return(out)
}

## The regression on which the plot is based 
reg02 <- lm(OUTsupport~TO*Tcond,data=subset(dbrazil,CCTbeneficiary==F&attention==T))
print(summary(reg02))

## Calculations for the plot
other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]

# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(dbrazil$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(dbrazil$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(d$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA

# Figure 1(A) in the paper
pdf(file=paste(savedir,"/fig-Brazil-P01S02-3x2barplot.pdf",sep=""))
par(mar=c(5,4.5,2,1))
tmp <- barplot(to.plot,beside=T,ylim=c(1,7), xpd=F)
polygon(x=c(3.5,6.5,6.5,3.5),y=c(0,0,8,8),border=NA,col=gray(0.8))
mtext(text="Support for Social Transfer Program",side=2,line=3)
#abline(h=to.plot[,1],lty=1,col=gray(0.8))
par(new=T)
barplot(to.plot,beside=T,
		ylim=if(max(predict(reg02))<1){c(0,1)}else{c(1,7)}, 
		xpd=F,legend.text= c("Non-conditional","Conditional"),
		args.legend = list(bty="n",cex=1.2))	
segments(x0=tmp[2,],x1=tmp[2,],
		y0=Pcond$lwr,
		y1=Pcond$upr,col=1,lwd=1)
axis(side=1,at=c(tmp[,c(2,1,3)]),
		tick=F,
		line=-1,
		labels=paste("N=", c(table(reg02$model$Tcond,reg02$model$TO)),sep=""))
text(x=tmp[2,],y=Pcond[,"upr"]
	,labels=paste("Effect = ",sprintf("%01.2f",effectssd),
			"sd\n(p=",
			sprintf("%01.2f",ps),")",sep="")
	,pos=3)
dev.off()





### This analysis is the same, but includes non-attentive respondents. 
### A variation of this figure is in the supplemental information ###

### Interaction Effects ####
robse <- function(x,sig=0.95){
	require(sandwich)
	robSE <- sqrt(diag(vcovHC(x,type="HC2")))
	terms <- length(coef(x))
	robP <- rep(NA,terms)
	for(jj in 1:terms){
			robP[jj] <- linearHypothesis(x,diag(terms)[jj,],vcov=vcovHC(x,type="HC2"))$Pr[2]}
	upr <- coef(x) + qnorm(sig) * robSE
	lwr <- coef(x) - qnorm(sig) * robSE
	out <- cbind(robse=robSE,upr=upr,lwr=lwr,robpval=robP)
	return(out)
}


reg02 <- lm(OUTsupport~TO*Tcond,data=subset(dbrazil,CCTbeneficiary==F))
print(summary(reg02))

other.order <- c("Racial","Control","Regional")
cond.order <- c("Unconditional","Conditional")

Puncond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("_Uncond",3),
			TO=other.order),
			interval="confidence"))
Pcond <- data.frame(predict(reg02,newdata=data.frame(
			Tcond=rep("Cond",3),
			TO=other.order),
			interval="confidence"))

to.plot <- rbind(Conditional=Pcond$fit,
				Unconditional=Puncond$fit)
colnames(to.plot)<-other.order
to.plot <- to.plot[cond.order,]

# standardized effect
effectssd <- c(racial=round((coef(reg02)[4]+coef(reg02)[5])/sd(dbrazil$OUTsupport,na.rm=T),2),
			control=round((coef(reg02)[4])/sd(dbrazil$OUTsupport,na.rm=T),2),
			regional=round((coef(reg02)[4]+coef(reg02)[6])/sd(d$OUTsupport,na.rm=T),2) )#effects in SD
	
# p-value of conditional relative unconditional		
ps <- round(c(racial=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Racial"))$p.value,
		control=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Control"))$p.value,
		regional=t.test(OUTsupport~Tcond,data=subset(reg02$model,TO=="Regional"))$p.value),3)
 
# p-value of premium in control relative to racial and regional
pps <-  round(robse(reg02)[c(5,4,6),"robpval"],3)
pps[2] <- NA

pdf(file=paste(savedir,"/fig-Brazil-P01S02SI-3x2barplot.pdf",sep=""))
par(mar=c(5,4.5,2,1))
tmp <- barplot(to.plot,beside=T,ylim=c(1,7), xpd=F)
polygon(x=c(3.5,6.5,6.5,3.5),y=c(0,0,8,8),border=NA,col=gray(0.8))
mtext(text="Support for Social Transfer Program",side=2,line=3)
#abline(h=to.plot[,1],lty=1,col=gray(0.8))
par(new=T)
barplot(to.plot,beside=T,
		ylim=if(max(predict(reg02))<1){c(0,1)}else{c(1,7)}, 
		xpd=F,legend.text= c("Non-conditional","Conditional"),
		args.legend = list(bty="n",cex=1.2))	
segments(x0=tmp[2,],x1=tmp[2,],
		y0=Pcond$lwr,
		y1=Pcond$upr,col=1,lwd=1)
axis(side=1,at=c(tmp[,c(2,1,3)]),
		tick=F,
		line=-1,
		labels=paste("N=", c(table(reg02$model$Tcond,reg02$model$TO)),sep=""))
text(x=tmp[2,],y=Pcond[,"upr"]
	,labels=paste("Effect = ",sprintf("%01.2f",effectssd),
			"sd\n(p=",
			sprintf("%01.2f",ps),")",sep="")
	,pos=3)
dev.off()
